function [Ntop] = TEC_TOP(h, NmF2, hmF2, H0)
%
% Topside Electron Density
%
%DESCRIPTION:
%This function computes the Topside Electron Density.
%
%PROTOTYPE:
% [Ntop] = TEC_TOP(h, NmF2, hmF2, H0)
%
%--------------------------------------------------------------------------
% INPUTS:
%   h          [1x1]       Height                    [km]
%   NmF2       [1x1]       F2-Layer Max. Density     [1e11 m-3]
%   hmF2       [1x1]       F2-Layer Max. Den. Height [km]
%   H0         [1x1]       Scale Height              [km]
%--------------------------------------------------------------------------
% OUTPUTS:
%   Ntop       [1x1]       Topside Electron Dens.    [m-3]
%--------------------------------------------------------------------------
%
%NOTES:
% (none)
%
%CALLED FUNCTIONS:
% (none)
%
%UPDATES:
% (none)
%
%REFERENCES:
% [1] "Ionospheric Correction Algorithm for Galileo Single-Frequency Users"
%      - European GNSS (Galileo) Open Service
% [2] "Electron Density Models and Data for Transionospheric Radio
%      Propagation" - Report ITU-R P.2297-1 (05/2019)
%
%AUTHOR(s):
%Luigi De Maria, Matteo D'Addazio, 2022
%

%% Main Code

%Constants
g = 0.125;
r = 100;

%Arguments
Dh = h - hmF2;
z = (Dh) / (H0 * (1 + (r*g*Dh) / (r*H0 + h*Dh)));

%Exponential
ea = exp(z);

%Electron Density
if ea > 1e11        %[2] does not report such condition and the 1e11 terms
    Ntop = 4 * NmF2/ea * 1e11;
else
    Ntop = 4 * NmF2 * ea/(1+ea)^2 * 1e11;
end

end